/* Copyright (C) 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012 Mario Orsi
   This file is part of Brahms.
   Brahms is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
   Brahms is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
   You should have received a copy of the GNU General Public License along with Brahms.  If not, see <http://www.gnu.org/licenses/>. */

/* topology.c -- loading the constant attributes of each site of the system in dedicated data structures */

#include "dataStructs.h"

extern Site *site;  /* array of Site, represents the state of the system */
extern const VecR region;
extern real **sig, **eps, **sigSquared, **productElectroMoments, *mass, *inertia;
extern real sigWat, epsWat;
extern const real extTemperature, deltaT;
extern const real sigChol, sigAmine, sigPhos, sigGly, sigEst, sigTail, sigCharge; 
extern const real epsChol, epsAmine, epsPhos, epsEst, epsGly, epsTail, epsCharge;
extern const real massChol, massAmine, massPhos, massGly, soluteMass, massEst, massTail, massWat;
extern const real inertiaWat, inertiaEst, inertiaGly;
extern const real cholCharge, phosCharge, glyDipole, estDipole, watDipole; 
extern const int nSites, nTypes, nLipids, nSolutes, zConstraint, nWaters;

real **CreateRealMatrix( int, int ); /* external function declaration */

// global variables
const real hBondWatEst, hBondWatPhos, hBondWatGly, hBondWatWat;
const real hBondAmineWat, hBondAminePhos, hBondAmineEst, hBondAmineGly; 
LipidStruct dopcStruct, dopeStruct;
real *electroMoment, *sigma, *epsilon, *charge, *dipole;
real **productElectroMoments; /* product of electrostatic moments, i.e. charge*charge charge*dipole or dipole*dipole */
VecR solInertia; // principal moments of inertia of the solute (rigid) molecule

void SetLipidTopologies() {
  dopcStruct.nBeads = 15;
  dopeStruct.nBeads = 15;
}

void LoadTopology( FILE *logFilePtr )
{
  int type, i, j, counter = 1;
  real rTemperature, rDensity, rDipole, rTimeStep;
  FILE *fp;
  
  fprintf(logFilePtr, "LoadTopology... "); fflush(logFilePtr);

  fp = fopen( "topology.log", "a" );
  
  sig = CreateRealMatrix( nTypes, nTypes ); 
  sigSquared = CreateRealMatrix( nTypes, nTypes ); 
  eps = CreateRealMatrix( nTypes, nTypes ); 
  productElectroMoments = CreateRealMatrix( nTypes, nTypes ); 

  AllocMem( sigma, nTypes, real ); /* allocating arrays */
  AllocMem( epsilon, nTypes, real );
  AllocMem( electroMoment, nTypes, real );
  AllocMem( charge, nTypes, real );
  AllocMem( dipole, nTypes, real );

  for ( type = 0; type < nTypes; type++ ) { /* fill in arrays */
    switch ( type ) {
    case WATER_TYPE: /* SSD water */ 
      sigma[ type ] = sigWat;
      epsilon[ type ] = epsWat;
      electroMoment[ type ] = watDipole;
      charge[ type ] = 0;
      dipole[ type ] = watDipole;
      break;
    case CHOLINE_TYPE: /* choline */
      sigma[ type ] = sigChol;
      epsilon[ type ] = epsChol;
      electroMoment[ type ] = cholCharge;
      charge[ type ] = cholCharge;
      dipole[ type ] = 0;
      break;
    case AMINE_TYPE: /* amine */
      sigma[ type ] = sigAmine;
      epsilon[ type ] = epsAmine;
      electroMoment[ type ] = cholCharge; /* amine (in DOPE) has the same charge as choline (in DOPC) */
      charge[ type ] = cholCharge;
      dipole[ type ] = 0;
      break;
    case PHOSPHATE_TYPE: /* phosphate */
      sigma[ type ] = sigPhos;
      epsilon[ type ] = epsPhos;
      electroMoment[ type ] = phosCharge;
      charge[ type ] = phosCharge;
      dipole[ type ] = 0;
      break;
    case GLYCEROL_TYPE: 
      sigma[ type ] = sigGly;
      epsilon[ type ] = epsGly;
      electroMoment[ type ] = glyDipole;
      charge[ type ] = 0;
      dipole[ type ] = glyDipole;
      break;
    case ESTER_TYPE: 
      sigma[ type ] = sigEst;
      epsilon[ type ] = epsEst;
      electroMoment[ type ] = estDipole;
      charge[ type ] = 0;
      dipole[ type ] = estDipole;
      break;
    case TAIL_TYPE: 
      sigma[ type ] = sigTail;
      epsilon[ type ] = epsTail;
      electroMoment[ type ] = 0;
      charge[ type ] = 0;
      dipole[ type ] = 0;
      break;
    case CHARGE_TYPE: 
      sigma[ type ] = sigCharge;
      epsilon[ type ] = epsCharge;
      electroMoment[ type ] = 0;
      charge[ type ] = 0;
      dipole[ type ] = 0;
      break;
    default: if ( !zConstraint ) printf( "\nERROR in LoadTopolgy() - got to the switch default case!\n\n" );
    }
  }

  /* compile the interaction matrices according to Lorentz-Berthelot mixing rules - ref: AllenTildesley, p.21 */
  for ( i = 0; i < nTypes; i++ ) {
    for ( j = 0; j < nTypes; j++ ) {
      sig[ i ][ j ] = 0.5 * ( sigma[ i ] + sigma[ j ] );
      sigSquared[ i ][ j ] = Sqr( sig[ i ][ j ] );
      eps[ i ][ j ] = sqrt( epsilon[ i ] * epsilon[ j ] );
      /* electrostatics pair-magnitudes */
      productElectroMoments[ i ][ j ] = electroMoment[ i ] * electroMoment[ j ];
      fprintf( fp, "item %6d: sig[ %d ][ %d ] = %.4f, sigSquared[ %d ][ %d ] = %.4f, eps[ %d ][ %d ] = %.4f \n",
	       counter, i, j, sig[ i ][ j ], i, j, sigSquared[ i ][ j ], i, j, eps[ i ][ j ] );
      fprintf( fp, "item %6d: productElectroMoments[ %d ][ %d ] = %.4f, \n", counter, i, j, productElectroMoments[ i ][ j ] );
      ++counter;
    }
  }

  if ( nWaters ) {
    mass[ WATER_TYPE ] = massWat;
    inertia[ WATER_TYPE ] = inertiaWat; // from input
    fprintf( fp, "mass[ WATER_TYPE ] = %2.0f u\n", mass[ WATER_TYPE ] );
    fprintf( fp, "inertia[ WATER_TYPE ] = %.4f u * nm^2\n", inertia[ WATER_TYPE ] );
    /* WATER <--> WATER hydrogen-bonding  */
    eps[ WATER_TYPE ][ WATER_TYPE ] *= hBondWatWat; 
    epsWat *= hBondWatWat; 
  }

  if ( nLipids ) {
    /* loading masses */
    mass[ CHOLINE_TYPE ] = massChol;
    mass[ AMINE_TYPE ] = massAmine;
    mass[ PHOSPHATE_TYPE ] = massPhos;
    mass[ GLYCEROL_TYPE ] = massGly;
    mass[ ESTER_TYPE ] = massEst;
    mass[ TAIL_TYPE ] = massTail;
    /* computing principal moments of inertia */
    /*     inertia[ GLYCEROL_TYPE ] = massGly * Sqr(sigGly) / 10; */
    /*     inertia[ ESTER_TYPE ]    = massEst * Sqr(sigEst) / 10; */
    inertia[ GLYCEROL_TYPE ] = inertiaGly;
    inertia[ ESTER_TYPE ]    = inertiaEst;
    fprintf( fp, "inertia[ GLYCEROL_TYPE ] = %.4f u * nm^2\n", inertia[ GLYCEROL_TYPE ] );
    fprintf( fp, "inertia[ ESTER_TYPE ]    = %.4f u * nm^2\n", inertia[ ESTER_TYPE ] );

    /* WATER --> PHOSPHATE hydrogen-bonding  */
    eps[ WATER_TYPE ][ PHOSPHATE_TYPE ] *= hBondWatPhos; 
    eps[ PHOSPHATE_TYPE ][ WATER_TYPE ] *= hBondWatPhos; 

    /* WATER <-- AMINE hydrogen-bonding  */
    eps[ WATER_TYPE ][ AMINE_TYPE ] *= hBondAmineWat; 
    eps[ AMINE_TYPE ][ WATER_TYPE ] *= hBondAmineWat; 

    /* WATER --> GLYCEROL hydrogen-bonding  */
    eps[ WATER_TYPE ][ GLYCEROL_TYPE ]  *= hBondWatGly;
    eps[ GLYCEROL_TYPE ][ WATER_TYPE ]  *= hBondWatGly;

    /* WATER --> ESTER hydrogen-bonding  */
    eps[ WATER_TYPE ][ ESTER_TYPE ]  *= hBondWatEst; 
    eps[ ESTER_TYPE ][ WATER_TYPE ]  *= hBondWatEst; 
    
    /* PHOSPHATE <-- AMINE hydrogen-bonding  */
    eps[ AMINE_TYPE ][ PHOSPHATE_TYPE ] *= hBondAminePhos; 
    eps[ PHOSPHATE_TYPE ][ AMINE_TYPE ] *= hBondAminePhos; 

    /* GLYCEROL <-- AMINE hydrogen-bonding  */
    eps[ AMINE_TYPE ][ GLYCEROL_TYPE ]  *= hBondAmineGly; 
    eps[ GLYCEROL_TYPE ][ AMINE_TYPE ]  *= hBondAmineGly; 

    /* ESTER <-- AMINE hydrogen-bonding  */
    eps[ AMINE_TYPE ][ ESTER_TYPE ]  *= hBondAmineEst; 
    eps[ ESTER_TYPE ][ AMINE_TYPE ]  *= hBondAmineEst; 
  }

  // just checking 
  for ( i = 0; i < nTypes; i++ ) {
    for ( j = 0; j < nTypes; j++ ) {
      counter = 1;
      fprintf( fp, "energy parameter %6d: eps[ %d ][ %d ] = %.4f kJ/mol\n", counter, i, j, eps[ i ][ j ]);
      ++counter;
    }
  }
    
  if (nSolutes) mass[ SOLUTE_TYPE ] = soluteMass;
    
  if ( zConstraint || nSolutes ) {
    printf("*** loading solute mass into the mass[] array...\n");
    mass[ SOLUTE ] = soluteMass;
    printf("*** mass[SOLUTE] = %7.3f\n", mass[ SOLUTE ] );
  }

  fprintf( fp, "\nReduced values for water (Stockmayer fluid):\n" );
  // From Appendix B of 'Computer simulation of liquids', Allen & Tildesley, 1987
  rTemperature = kB_IN_BRAHMS_UNITS * extTemperature / (hBondWatWat*epsWat);
  rDensity = (nWaters / (region.x*region.y*region.z)) * Cube( sigWat ); 
  rDipole = watDipole / sqrt( 4 * PI * EPSILON0_IN_BRAHMS_UNITS * (hBondWatWat*epsWat) * Cube( sigWat ) );
  rTimeStep = sqrt( epsWat / ( mass[ WATER_TYPE ] * Sqr( sigWat ) ) ) * deltaT;
  fprintf( fp, "Temperature* = %.2f\n", rTemperature );
  fprintf( fp, "Density* = %.2f\n", rDensity );
  fprintf( fp, "Dipole* = %.2f\n", rDipole );
  fprintf( fp, "TimeStep* = %.4f\n\n", rTimeStep );

  
      
  fclose( fp );
  fprintf(logFilePtr, "done\n"); fflush(logFilePtr);
}

